Tmap

We will map the HFMD cases notified in Pulau Pinang in 2022

Packages installation and load

library(tmap) ## thematic maps in R
library(sf) ## simple feature creation for mapping

Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(readxl) ## to load data
library(here) ## read data

here() starts at C:/Users/ACER/Desktop/Rconf_Tmap

library(dplyr) ##data wrangling

Warning: package 'dplyr' was built under R version 4.2.3

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

Part 1

Importing the data

data("World") #open source map in SF package

Now that we have the map, we can plot the “World” map by using function tm_shape. tm_shape will determine the layer of map that we wan to plot as base map.

tm_polygon determine which variable that we want to show on the chloropleth map

tm_shape(World) +
  tm_polygons("life_exp")

World map produced, however the breaks and scales used were not really show any differences

We need to modify it with (break = ) functions in the package.

tm_shape(World) +
  tm_polygons("life_exp",
              breaks = c (0, 30, 40, 50, 60,70,80, Inf))

Lets customize the map to make it more informative and reader’s friendly

tm_shape(World) +
  tm_polygons("life_exp",
              title = "Life Expectancy" , ## add title to the legend
              breaks = c (0, 30, 40, 50, 60,70,80, Inf)) +
    tm_borders() ## add borders on the map 

Warning: One tm layer group has duplicated layer types, which are omitted. To
draw multiple layers of the same type, use multiple layer groups (i.e. specify
tm_shape prior to each of them).

We can view it in dynamic map view Just need to activate and shif the plotting mode into view mode

Activate “view” mode

tmap_mode("view")

tmap mode set to interactive viewing

Rerun the map We can click on the map or zoom-in the area/ country in this example, population density information will be appear once we click any country

Plot the map

tm_shape(World) +
  tm_polygons("life_exp",
              title = "Life Expectancy" , ## add title to the legend
              breaks = c (0, 30, 40, 50, 60,70,80, Inf)) +
    tm_borders() ## add borders on the map 

Warning: One tm layer group has duplicated layer types, which are omitted. To
draw multiple layers of the same type, use multiple layer groups (i.e. specify
tm_shape prior to each of them).

we can edit the pop-up to make it more readers-friendly

tm_shape(World) +
  tm_polygons("life_exp",
              title = "Life Expectancy" , ## add title to the legend
               breaks = c (0, 30, 40, 50, 60,70,80, Inf),
              id = "name",
              popup.vars = c("Country" = "name","Life_Expectancy" = "life_exp", "Pop_Density" = "pop_est_dens")) +
    tm_borders() ## add borders on the map 

Warning: One tm layer group has duplicated layer types, which are omitted. To
draw multiple layers of the same type, use multiple layer groups (i.e. specify
tm_shape prior to each of them).

Part 2

Lets map the disease

We need base map (area of interest) and data of the cases including the coordinates

Load the case data

Case <- read_excel("CASE.xlsx")

glimpse(Case)

Rows: 5,924
Columns: 9
$ WEEK        <dbl> 20, 20, 31, 24, 23, 17, 19, 15, 10, 19, 4, 26, 20, 20, 19,…
$ AGE         <dbl> 2, 5, 8, 5, 2, 4, 4, 1, 1, 3, 4, 6, 4, 5, 3, 3, 11, 2, 7, …
$ ADDRESS     <chr> "962 A MUKIM 20 PENANTI", "TBP 71P, TANAH LIAT, 14000, BM"…
$ RACE        <chr> "Melayu", "Melayu", "Melayu", "Melayu", "Melayu", "Melayu"…
$ NATIONALITY <chr> "Warganegara", "Warganegara", "Warganegara", "Warganegara"…
$ GENDER      <chr> "Lelaki", "Perempuan", "Lelaki", "Perempuan", "Perempuan",…
$ DAERAH      <chr> "SEBERG PERAI TENGAH", "SEBERG PERAI TENGAH", "SEBERG PERA…
$ LATITUDE    <chr> "5.406580000000076", "5.3805049999526275", "5.397249988930…
$ LONGITUDE   <chr> "100.47596000000004", "100.46462502352176", "100.403379965…

Convert the case data to an SF object

Case_1 <- st_as_sf(Case,coords = c("LONGITUDE","LATITUDE"), 
                    crs = 4326)  ## give coordinate referral system 
Case_1 <- st_transform(Case_1, 3168) ## Transform the CRS into Kertau format

Plot the case

tm_shape(Case_1) + 
  tm_dots(col = "black")

Load the base map

Penang <- st_read(here("Penang", 
                  "mukim2000_penang1.shp"))

Reading layer `mukim2000_penang1' from data source 
  `C:\Users\ACER\Desktop\Rconf_Tmap\Penang\mukim2000_penang1.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 83 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 243501.5 ymin: 567201.1 xmax: 285076.2 ymax: 618503.2
CRS:           NA

Give CRS to base map (similar to case data)

st_crs(Penang)<-3168  ## Transform the CRS into Kertau format

Plotting mode

tmap_mode("plot") ## plotting mode

tmap mode set to plotting

Now that we have both cases and base map, we can plot the cases on the map

tm_shape(Penang) +  ## The base map layer
  tm_polygons("NAMA_DP", 
              title = "HFMD Cases in Pulau Pinang 2022",
              title.size = 10) +
  tm_shape(Case_1) +  ## the case layer
  tm_dots(col = "black")+ ## use tm_dots to represent each case
  tm_layout(legend.outside = TRUE)

We can plot the cases according to specific district

We just need to adjust the layer

tm_shape(Penang[Penang$NAMA_DP == "SP.TENGAH",])+ ## specify the district
  tm_polygons("NAMMUK",
              title = "HFMD Cases in SP.TENGAH in 2022")+
  tm_shape(Case_1[Case_1$DAERAH == "SEBERG PERAI TENGAH",])+ ## specify the district in 2nd ayer
  tm_dots (col = "black") +  ## use tm_dots to represent each case
  tm_layout(legend.outside = TRUE)

We also can make direct comparison between groups (eg: gender) by using tm_facet function

Tm_Facet function

tm_shape(Penang) +  ## The base map layer
  tm_polygons("NAMA_DP", 
              title = "HFMD Cases in Pulau Pinang 2022",
              title.size = 10) +
  tm_shape(Case_1) +  ## the case layer
  tm_dots(col = "black")+   ## use tm_dots to represent each case
  tm_layout(legend.outside = TRUE) +
  tm_facets("GENDER", ncol = 2) 

Now, we can try it for areal data too (mapping the incidence rate)

Lets load the data

Case_2 <- read_excel( "Case_Sub.xlsx")
Case_3 <- merge(Penang, Case_2, by = c("NAMA_DP", "NAMMUK"))  
Incidence <- Case_3 %>% mutate(Incidence_Rate = (n/PENDUDUK)*1000)

glimpse(Incidence)

Rows: 80
Columns: 7
$ NAMA_DP        <chr> "BARAT DAYA", "BARAT DAYA", "BARAT DAYA", "BARAT DAYA",…
$ NAMMUK         <chr> "MUKIM 1(PANTAI ACHEH)", "MUKIM 10(BKT.RELAU)", "MUKIM …
$ NAMA_NG        <chr> "PULAU PINANG", "PULAU PINANG", "PULAU PINANG", "PULAU …
$ PENDUDUK       <dbl> 4422, 2195, 11734, 94740, 3234, 2035, 1822, 145, 6261, …
$ n              <dbl> 17, 110, 172, 606, 23, 11, 21, 8, 38, 9, 145, 10, 23, 2…
$ geometry       <MULTIPOLYGON [m]> MULTIPOLYGON (((247469.7 60..., MULTIPOLYG…
$ Incidence_Rate <dbl> 3.844414, 50.113895, 14.658258, 6.396453, 7.111936, 5.4…

Plot the incidence

tm_shape(Incidence) +
  tm_polygons("Incidence_Rate",
              title = "Incidence of HFMD Pulau Pinang 2022",
              title.size = 1,
             breaks = c (0,6,9,12,18, Inf)) +
    tm_compass(type = "4star", size = 2, position = c("left", "top")) + 
  tm_layout(legend.outside = TRUE)

Lets make it dynamic

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(Incidence) +
  tm_polygons("Incidence_Rate",
              title = "Incidence of HFMD Pulau Pinang 2022",
              title.size = 1,
             breaks = c (0,6,9,12,18, Inf),
             id = "NAMA_DP",
              popup.vars = c("Subdistrict" = "NAMMUK","Case" = "n", "Incidence" = "Incidence_Rate")) +
    tm_borders() ## add borders on the map 

Warning: One tm layer group has duplicated layer types, which are omitted. To
draw multiple layers of the same type, use multiple layer groups (i.e. specify
tm_shape prior to each of them).